rm(list=ls())
library(haven)
Warning: package ‘haven’ was built under R version 4.3.3
# lettura dati
setwd("C:/Users/dario/Documents/Projects/Master/Geospatial/")
stat <- read_dta("data/slave_trade_QJE.dta")
# rinonimo la colonna del PIL Log per capita GDP - from Maddison (2003)
colnames(stat)[colnames(stat) == "ln_maddison_pcgdp2000"] = "ln_pcgdp"
colnames(stat)[colnames(stat) == "ln_coastline_area"] = "ln_coast_area"
colnames(stat)[colnames(stat) == "ln_avg_gold_pop"] = "ln_gold_pop"
colnames(stat)[colnames(stat) == "ln_avg_oil_pop"] = "ln_oil_pop"
colnames(stat)[colnames(stat) == "ln_avg_all_diamonds_pop"] = "ln_diamonds_pop"
colnames(stat)[colnames(stat) == "atlantic_distance_minimum"] = "atlantic_dist_min"
colnames(stat)[colnames(stat) == "indian_distance_minimum"] = "indian_dist_min"
colnames(stat)[colnames(stat) == "saharan_distance_minimum"] = "saharan_dist_min"
colnames(stat)[colnames(stat) == "red_sea_distance_minimum"] = "red_sea_dist_min"
# rimuovo variabili non necessarie
delenda = c(#'abs_latitude', # Absolute latitude
#'longitude', #Longitude
'ln_export_pop', #Log total slave exports normalized by historic population
'island_dum', #Indicator variable for small islands
#'region_n', # Region indicator: North
'region_s', # Region indicator: South
'region_w', # Region indicator: West
'region_e', # Region indicator: East
'region_c', # Region indicator: Central
'ln_pop_dens_1400', # Log population density in 1400
'ethnic_fractionalization', # Ethnic fractionalization
'state_dev', # State development
'land_area') # Land area in millions of square kms
stat[,delenda] = NULL
library(sf)
Warning: package ‘sf’ was built under R version 4.3.3Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
# lettura shapefile
world = st_read('data/countries/50m/ne_50m_admin_0_countries.shp')
Reading layer `ne_50m_admin_0_countries' from data source
`C:\Users\dario\Documents\Projects\Master\Geospatial\data\countries\50m\ne_50m_admin_0_countries.shp'
using driver `ESRI Shapefile'
Simple feature collection with 242 features and 168 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
Geodetic CRS: WGS 84
africa = world[which(world$CONTINENT=='Africa'),c('ISO_A3', 'NAME')]
par(mar=c(0,0,0,9))
coords <- st_coordinates(st_centroid(africa))
Warning: st_centroid assumes attributes are constant over geometries
plot(st_geometry(africa), border='gray', col='snow')
text(coords[, "X"], coords[, "Y"], labels=africa$ISO_A3, cex=0.5)
legend(x=st_bbox(africa)$xmax+1, y=st_bbox(africa)$ymax,
legend = paste(africa$ISO_A3, ' ', africa$NAME), cex=0.5, ncol=2, xpd=T)
par(mar=c(5, 4, 4, 2)+0.1)
Esiste un legame causale tra le deportazioni degli schiavi
(ln_export_area) ed il PIL procapite
(ln_pcgdp) degli odierni Stati africani?
par(mar=c(4.5,4,0.5,0.5))
plot(x=stat$ln_export_area, y=stat$ln_pcgdp, pch=16, bty='n', cex=0.7,
xlab='ln_export_area', ylab='ln_pcgdp', col='firebrick')
text(x=stat$ln_export_area, y=stat$ln_pcgdp+0.1, labels=stat$isocode, cex=0.5)
grid()
par(mar=c(5,4,4,2)+0.1)
#rimuovo le isole
delenda = c('SYC', # Seychelles
'STP', # Sao Tome & Principe
'MUS', # Mauritius
'MDG', # Madagascar
'CPV', # Cape Verde Islands
'COM' # Comoros
)
stat = stat[!(stat$isocode %in% delenda), ]
africa = africa[!(africa$ISO_A3 %in% delenda),]
Rispetto a quando i dati sul PIL sono stati raccolti (fine anni ’90), la conformazione politica dell’Africa è cambiata; alcuni Stati si sono separati, per cui è necessario riunificare le loro mappe attuali:
merge_countries <- function(sf_data, iso1, iso2, merged_iso, merged_name)
{
country1 = sf_data[which(sf_data$ISO_A3 == iso1), ]
country2 = sf_data[which(sf_data$ISO_A3 == iso2), ]
merged_geometry = st_union(st_geometry(country1), st_geometry(country2))
merged_data = data.frame(ISO_A3=merged_iso, NAME=merged_name, geometry=merged_geometry)
merged_sf = st_sf(merged_data, crs=st_crs(sf_data))
sf_data_filtered = sf_data[!(sf_data$ISO_A3 %in% c(iso1, iso2)), ]
sf_data_merged = rbind(sf_data_filtered, merged_sf)
return(sf_data_merged)
}
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])
(africa[tofix_shape,c('NAME', 'ISO_A3')])
Simple feature collection with 5 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -17.09878 ymin: -13.45381 xmax: 48.93857 ymax: 27.65645
Geodetic CRS: WGS 84
NAME ISO_A3 geometry
52 S. Sudan SSD MULTIPOLYGON (((33.97607 4....
59 Somaliland -99 MULTIPOLYGON (((48.93857 11...
105 W. Sahara ESH MULTIPOLYGON (((-8.817773 2...
174 Eritrea ERI MULTIPOLYGON (((36.52432 14...
192 Dem. Rep. Congo COD MULTIPOLYGON (((30.75117 -8...
# Rinonimo Dem. Rep. Congo (COD) in Democratic Republic of Congo (ZAR)
africa$ISO_A3[which(africa$ISO_A3=='COD')] = 'ZAR'
africa$NAME[which(africa$ISO_A3=='COD')] = 'Democratic Republic of Congo'
# Eritrea -> Etiopia
africa = merge_countries(africa, 'ERI', 'ETH', 'ETH', 'Ethiope')
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])
(africa[tofix_shape, c('NAME', 'ISO_A3')])
Simple feature collection with 3 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -17.09878 ymin: 3.490723 xmax: 48.93857 ymax: 27.65645
Geodetic CRS: WGS 84
NAME ISO_A3 geometry
52 S. Sudan SSD MULTIPOLYGON (((33.97607 4....
59 Somaliland -99 MULTIPOLYGON (((48.93857 11...
105 W. Sahara ESH MULTIPOLYGON (((-8.817773 2...
# W. Sahara -> Morocco
africa = merge_countries(africa, 'ESH', 'MAR', 'MAR', 'Morocco')
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])
(africa[tofix_shape, c('NAME', 'ISO_A3')])
Simple feature collection with 2 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 24.14736 ymin: 3.490723 xmax: 48.93857 ymax: 12.2231
Geodetic CRS: WGS 84
NAME ISO_A3 geometry
52 S. Sudan SSD MULTIPOLYGON (((33.97607 4....
59 Somaliland -99 MULTIPOLYGON (((48.93857 11...
# Somaliland -> Somalia
africa = merge_countries(africa, '-99', 'SOM', 'SOM', 'Somalia')
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])
(africa[tofix_shape, c('NAME', 'ISO_A3')])
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 24.14736 ymin: 3.490723 xmax: 35.26836 ymax: 12.2231
Geodetic CRS: WGS 84
NAME ISO_A3 geometry
52 S. Sudan SSD MULTIPOLYGON (((33.97607 4....
# S. Sudan -> Sudan
africa = merge_countries(africa, 'SSD', 'SDN', 'SDN', 'Sudan')
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])
(africa[tofix_shape,c('NAME', 'ISO_A3')])
Simple feature collection with 0 features and 2 fields
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84
[1] NAME ISO_A3 geometry
<0 rows> (or 0-length row.names)
Rimuoviamo le isole appartenenti ad alcuni Stati (in particolare un arcipelago molto a Sud sotto la nazionalità del Sudafrica)
### Rimozione isole
remove_islands <- function(sf, iso_a3)
{
country = sf[which(sf$ISO_A3 == iso_a3), ]
if (nrow(country) > 0)
{
geom_type <- st_geometry_type(country)[1] # Get the geometry type of the FIRST feature
if (geom_type %in% c("MULTIPOLYGON"))
{
#country_polygons <- st_cast(country, "POLYGON")
country_polygons <- suppressWarnings(st_cast(country, "POLYGON")) # Suppress warning here
areas = st_area(country_polygons)
largest_polygon = country_polygons[which.max(areas),]
# Extract the geometry of the largest polygon
country_geom = st_geometry(largest_polygon)[[1]]
# Create a new sf object for the mainland
country_mainland_sf = st_sf(data.frame(ISO_A3 = iso_a3, NAME = country$NAME[1]),
geometry = st_sfc(country_geom, crs = st_crs(country)))
sf_filtered = sf[!(sf$ISO_A3 %in% iso_a3), ]
sf_final <- rbind(sf_filtered, country_mainland_sf)
return(sf_final)
}
}
return(sf)
}
# Rimozione delle isole di ogni stato
delenda = c('TUN', # Tunisia
'TZA', # Tanzania
'ZAF', # S. Africa
'SLE', # Sierra Leone
'MRT', # Mauritania
'MWI', # Malawi
'GNB', # Guinea-Bissau
'GNQ', # Eq. Guinea
'ETH') # Ethiope
par(mfrow = c(1, 2), mar=c(0,0,1,0))
for (c in delenda)
{
plot(st_geometry(africa[which(africa$ISO_A3==c),'ISO_A3']), main=c, col='snow')
africa = remove_islands(africa, c)
plot(st_geometry(africa[which(africa$ISO_A3==c),'ISO_A3']), main=c, col='snow')
}
par(mfrow = c(1,1), mar=c(5, 4, 4, 2)+0.1)
names(africa)[1:2] = names(stat[1:2])
stat$country = NULL #non mi serve più
ds = merge(africa, stat, by = "isocode")
coords <- st_coordinates(st_centroid(ds))
Warning: st_centroid assumes attributes are constant over geometries
par(mar=c(0,0,0,8))
plot(st_geometry(ds), border='gray', col='snow')
text(coords[, "X"], coords[, "Y"], labels=ds$isocode, cex=0.5)
legend(x=st_bbox(ds)$xmax+1, y=st_bbox(ds)$ymax,
legend = paste(ds$isocode, ' ', ds$country), cex=0.5, ncol=2, xpd=T)
par(mar=c(5, 4, 4, 2)+0.1)
#utm <- st_crs("+proj=utm +zone=33 +north +datum=WGS84 +units=m +no_defs")
utm <- st_crs(2312)
sf_utm <- st_transform(ds, crs = utm)
st_crs(sf_utm)
Coordinate Reference System:
User input: EPSG:2312
wkt:
PROJCRS["Garoua / UTM zone 33N",
BASEGEOGCRS["Garoua",
DATUM["Garoua",
ELLIPSOID["Clarke 1880 (RGS)",6378249.145,293.465,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4197]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Cameroon - Garoua area."],
BBOX[8.92,12.9,9.87,14.19]],
ID["EPSG",2312]]
par(mfrow = c(1, 2), mar=c(0,0,2,0))
plot(st_geometry(ds), col = "lightcyan", main = "Original Projection")
plot(st_geometry(sf_utm), col = "lightyellow", main = "EPSG:2312 Projection")
par(mfrow = c(1, 1), mar=c(5, 4, 4, 2)+0.1)
# Lo salvo come shapefile per aprirlo in Geoda
# Save as Shapefile
st_write(sf_utm, "./data/slaveTrade.shp", driver="ESRI Shapefile", delete_layer=T)
Warning: Field names abbreviated for ESRI Shapefile driver
Deleting layer `slaveTrade' using driver `ESRI Shapefile'
Writing layer `slaveTrade' to data source
`./data/slaveTrade.shp' using driver `ESRI Shapefile'
Writing 46 features with 29 fields and geometry type Unknown (any).
# Lo converto in SpatialPolygonDataFrame per ri-usare codice prof
slaveTrade = as_Spatial(sf_utm)
isocode Country isocodecountry Country nameln_pcgdp Log per capita GDP - from Maddison (2003)ln_export_area Log total slave exports normalized by
land areacolony0 Colonizer indicator: not colonizedcolony1 Colonizer indicator: Britancolony2 Colonizer indicator: Francecolony3 Colonizer indicator: Portugalcolony4 Colonizer indicator: Belgiumcolony5 Colonizer indicator: Spaincolony6 Colonizer indicator: UNcolony7 Colonizer indicator: Italyabs_latitude Absolute latitudelongitude Longituderain_min Min of monthly average rainfall (mm)humid_max Max of monthly afternoon avg humidity
(%)low_temp Min of avg monthly low temp (C)ln_coastline_area Log (coastline/land_area)islam Percent Islamiclegor_fr Legal origin indicator: Frenchlegor_uk Legal origin indicator: Britishregion_n Region indicator: North Africaln_avg_gold_pop Log gold production per capitaln_avg_oil_pop Log oil production per capitaln_avg_all_diamonds_pop Log diamond production per
capitaatlantic_distance_minimum Minimum Atlantic distance
(000s of kms)indian_distance_minimum Minimum Indian distance (000s
of kms)saharan_distance_minimum Minimum Saharan distance (000s
of kms)red_sea_distance_minimum Minimum Red Sea distance (000s
of kms)library(sp)
library(spdep)
Warning: package ‘spdep’ was built under R version 4.3.3Loading required package: spData
Warning: package ‘spData’ was built under R version 4.3.3To access larger datasets in this package, install the spDataLarge package
with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Attaching package: ‘spData’
The following objects are masked _by_ ‘.GlobalEnv’:
coords, world
library(spatialreg)
Warning: package ‘spatialreg’ was built under R version 4.3.3Loading required package: Matrix
Attaching package: ‘spatialreg’
The following objects are masked from ‘package:spdep’:
get.ClusterOption, get.coresOption, get.mcOption, get.VerboseOption,
get.ZeroPolicyOption, set.ClusterOption, set.coresOption, set.mcOption,
set.VerboseOption, set.ZeroPolicyOption
library(RColorBrewer)
library(fields) # For the colorbar
Warning: package ‘fields’ was built under R version 4.3.3Loading required package: spam
Warning: package ‘spam’ was built under R version 4.3.3Spam version 2.11-0 (2024-10-03) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
Attaching package: ‘spam’
The following object is masked from ‘package:Matrix’:
det
The following objects are masked from ‘package:base’:
backsolve, forwardsolve
Loading required package: viridisLite
Warning: package ‘viridisLite’ was built under R version 4.3.3
Try help(fields) to get started.
library(gridExtra)
Warning: package ‘gridExtra’ was built under R version 4.3.3
matlab.like.hot <- function(n)
{
#my_gradient = c('gray25', 'red', 'yellow', 'floralwhite')
my_gradient = c('gray25', 'red', 'yellow', 'cornsilk')
#my_gradient = c('gray25', 'red', 'yellow', 'lightgoldenrodyellow')
#my_gradient = c('gray25', 'red', 'yellow', 'lightyellow')
#my_gradient = c('gray25', 'red', 'yellow', 'lemonchiffon')
#my_gradient = c('gray25', 'red', 'yellow', 'snow')
#my_gradient = c('gray25', 'red', 'yellow', 'ivory')
return (colorRampPalette(my_gradient, space = "Lab")(n))
}
plotMapData <- function(shp, vname, min_v=NA, max_v=NA)
{
n_colors <- 16
# spectral_colors = brewer.pal(11, "Spectral")
# spectral_r_palette <- colorRampPalette(rev(spectral_colors)) # Reverse for _r
palette = matlab.like.hot(n_colors) #spectral_r_palette(n_colors)
idv = which(names(shp)==vname)
vals = shp[[idv]]
if (is.na(min_v)) (min_v = min(vals))
if (is.na(max_v)) (max_v = max(vals))
if (min_v == max_v)
{
(min_v = floor(min(vals)))
(max_v = ceiling(max(vals)))
}
norm_v = ((vals - min_v) / (max_v - min_v))
norm_v[norm_v < 0] = 0
norm_v[norm_v > 1] = 1
pcol = palette[round(norm_v * (n_colors-1)) + 1]
par(mar = c(0.0, 0.0, 2.0, 3))
plot(shp, col = pcol, lwd=0.5, main=vname, border='lightskyblue4')
# Colorbar
image.plot(legend.only = TRUE, zlim = c(min_v, max_v),
col = palette,
axis.args=list(cex.axis=0.8))
par(mar=c(5,4,4,2) + 0.1)# Reset to default layout
# spplot(shp, vname, main=vname, lwd=0.5,
# col = 'lightskyblue4',
# col.regions = palette,
# par.settings = list(axis.line = list(col = NA)))
}
showData <- function(shp, col_name, min_v=NA, max_v=NA)
{
#par(mfrow = c(1, 2), mar = c(3, 4, 1.0, 0.5))
par(mfrow = c(1, 2), mar = c(2, 4, 1.0, 0.5))
h = hist(shp[[col_name]],main=col_name,xlab='',
col='lightskyblue3', border='lightskyblue4')#, breaks='FD',)
grid(lty='solid', lwd=0.5, col='white', nx=NA, ny=NULL)
plotMapData(shp, col_name)
par(mfrow=c(1,1), mar=c(5,4,4,2)+0.1)
}
#--------------------------------------------------------------------
#ln_maddison_pcgdp2000
showData(slaveTrade,'ln_pcgdp')
#ln_export_area
showData(slaveTrade,'ln_export_area')
# rain_min
showData(slaveTrade,'rain_min')
# humid_max
showData(slaveTrade,'humid_max')
# low_temp
showData(slaveTrade,'low_temp')
# ln_coastline_area
showData(slaveTrade,'ln_coast_area')
# islam
showData(slaveTrade,'islam')
# ln_avg_gold_pop
showData(slaveTrade,'ln_gold_pop')
# ln_avg_oil_pop
showData(slaveTrade,'ln_oil_pop')
# ln_avg_all_diamonds_pop
showData(slaveTrade,'ln_diamonds_pop')
# atlantic_distance_minimum
showData(slaveTrade,'atlantic_dist_min')
# indian_distance_minimum
showData(slaveTrade,'indian_dist_min')
# saharan_distance_minimum
showData(slaveTrade,'saharan_dist_min')
# red_sea_distance_minimum
showData(slaveTrade,'red_sea_dist_min')
#..........................................................
colony_codes = c('--', 'GB', 'FR', 'PT', 'BE', 'ES', 'UN', 'IT')
#colony_names = names(slaveTrade)[5:12]
colony_names = c('colony0', 'colony1', 'colony2', 'colony3', 'colony4', 'colony5', 'colony6', 'colony7')
colony_id = max.col(stat[,colony_names])
slaveTrade$colony = as.factor(colony_codes[colony_id])
legor_codes = c('FR', 'GB', '--')
legor_names = c('legor_fr', 'legor_uk')
legor_id = max.col(stat[,legor_names])
slaveTrade$legor = legor_codes[legor_id]
slaveTrade$legor[rowSums(stat[,legor_names]) == 0] = '--'
slaveTrade$legor = as.factor(legor_codes[legor_id])
library(gridExtra)
#pL =
spplot(slaveTrade, 'colony', main='colonizers', lwd=0.5,
col.regions = brewer.pal(n=length(colony_codes), name='Pastel2'), par.settings = list(axis.line = list(col = NA)))
# pR =
spplot(slaveTrade, 'legor', main='legislative origin', lwd=0.5,
col.regions = brewer.pal(n=length(legor_codes), name='Pastel1'), par.settings = list(axis.line = list(col = NA)))
#grid.arrange(pL, pR, ncol = 2)
library(corrplot)
Warning: package ‘corrplot’ was built under R version 4.3.3corrplot 0.95 loaded
#(preds = names(slaveTrade)[c(3:4,13:18, 20:26)])
preds = c('ln_export_area', 'abs_latitude', 'longitude',
'rain_min', 'humid_max', 'low_temp', 'ln_coast_area',
'islam', 'ln_gold_pop', 'ln_oil_pop', 'ln_diamonds_pop',
'atlantic_dist_min', 'indian_dist_min',
'saharan_dist_min', 'red_sea_dist_min')
corpredictor = cor(slaveTrade@data[,preds], method="pearson")
par(mar=c(0,0,0,0))
corrplot(corpredictor, type = "upper", col=rev(colorRampPalette(brewer.pal(n=11, name='RdBu'))(100)), tl.col = "black", tl.srt = 45, tl.cex=0.8)
par(mar=c(5,4,4,2)+0.1)
(corpredictor['abs_latitude', 'low_temp'])
[1] -0.7768648
(corpredictor['longitude', c('atlantic_dist_min', 'indian_dist_min', 'red_sea_dist_min')])
atlantic_dist_min indian_dist_min red_sea_dist_min
0.8145862 -0.7918568 -0.7828318
(corpredictor['islam', 'saharan_dist_min'])
[1] -0.7350223
(corpredictor['atlantic_dist_min', 'red_sea_dist_min'])
[1] -0.8306802
slaveTrade.nb <- poly2nb(slaveTrade)
slaveTrade.lw <- nb2listw(slaveTrade.nb, style = "W")
slaveTrade.nb.lag = nblag(slaveTrade.nb, maxlag=5)
summary(slaveTrade.nb)
Neighbour list object:
Number of regions: 46
Number of nonzero links: 202
Percentage nonzero weights: 9.546314
Average number of links: 4.391304
Link number distribution:
1 2 3 4 5 6 7 8 9
2 8 7 7 8 8 2 3 1
2 least connected regions:
17 23 with 1 link
1 most connected region:
44 with 9 links
summary(slaveTrade.lw)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 46
Number of nonzero links: 202
Percentage nonzero weights: 9.546314
Average number of links: 4.391304
Link number distribution:
1 2 3 4 5 6 7 8 9
2 8 7 7 8 8 2 3 1
2 least connected regions:
17 23 with 1 link
1 most connected region:
44 with 9 links
Weights style: W
Weights constants summary:
cnb = coordinates(slaveTrade)
par(mfrow=c(1,1), mar=c(0,0,1,0))
plot(slaveTrade, border='gray', lwd=0.5, main='prox. lag 1')
plot(slaveTrade.nb.lag[[1]], cnb, add=T, col='red', pch=16, cex=0.5, lwd=0.5)
#
par(mfrow=c(2,2), mar=c(0,0,1,0))
#
plot(slaveTrade, border='gray', main='prox. lag 2', lwd=0.5)
plot(slaveTrade.nb.lag[[2]], cnb, add=T, col='red', pch=16, cex=0.3, lwd=0.5)
#
plot(slaveTrade, border='gray', main='prox. lag 3', lwd=0.5)
plot(slaveTrade.nb.lag[[3]], cnb, add=T, col='red', pch=16, cex=0.3, lwd=0.5)
#
plot(slaveTrade, border='gray', main='prox. lag 4', lwd=0.5)
plot(slaveTrade.nb.lag[[4]], cnb, add=T, col='red', pch=16, cex=0.3, lwd=0.5)
#
plot(slaveTrade, border='gray', main='prox. lag 5', lwd=0.5)
plot(slaveTrade.nb.lag[[5]], cnb, add=T, col='red', pch=16, cex=0.3, lwd=0.5)
par(mfrow=c(2,2), mar=c(5, 4, 4, 2)+0.1)
library(pgirmess)
moran_list = list()
#cvars = names(slaveTrade)[c(3,4,13:26)]
cvars = c("ln_pcgdp", "ln_export_area", "rain_min", "humid_max", "low_temp", "ln_coast_area",
"islam",
"ln_gold_pop","ln_oil_pop", "ln_diamonds_pop",
"atlantic_dist_min", "indian_dist_min", "saharan_dist_min", "red_sea_dist_min")
for (varname in cvars)
{
vals = slaveTrade[[varname]]
resM = moran.test(vals, slaveTrade.lw)
resG = geary.test(vals, slaveTrade.lw)
moran_list = rbind(moran_list, c(varname, resM$estimate[1], resM$p.value, resG$estimate[1], resG$p.value))
corD<-correlog(coordinates(slaveTrade), slaveTrade[[varname]], method="Moran")
par(mar=c(4, 3, 2, 0))
barplot(corD[,'coef'], names.arg=1:length(corD[,'coef']), yaxt = "n",
main=paste('correlog. ', varname), xlab='lag', ylim=c(-1,1),
col='lightskyblue3', border='lightskyblue4')
axis(2, # 2 indicates the left axis (y-axis)
at = seq(-1,1,0.2), # Positions of the ticks
labels = seq(-1,1,0.2))
#grid(lty='solid', lwd=0.5, col='white', nx=NA, ny=NULL)
abline(h = seq(-1,1,0.2), col = "white", lty = 'solid', lwd=0.5) # Add gridlines at y_ticks
par(mar=c(5, 4, 4, 2)+0.1)
print(varname)
}
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "ln_pcgdp"
[1] "ln_export_area"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "rain_min"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "humid_max"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "low_temp"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "ln_coast_area"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "islam"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "ln_gold_pop"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "ln_oil_pop"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "ln_diamonds_pop"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "atlantic_dist_min"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "indian_dist_min"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "saharan_dist_min"
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
[1] "red_sea_dist_min"
print(moran_list)
Moran I statistic
[1,] "ln_pcgdp" "0.423212218627437" "5.84693606250667e-06"
[2,] "ln_export_area" "0.305517058870406" "0.000707523106463913"
[3,] "rain_min" "0.164251212033616" "0.0270892478720973"
[4,] "humid_max" "0.273673184682141" "0.00146370187499353"
[5,] "low_temp" "0.569670215279815" "3.870356674188e-09"
[6,] "ln_coast_area" "0.242017671557775" "0.00510666852032284"
[7,] "islam" "0.705157163497059" "7.96725428967743e-13"
[8,] "ln_gold_pop" "0.174254358907525" "0.0284356869275728"
[9,] "ln_oil_pop" "0.467275141106285" "7.95651993408283e-07"
[10,] "ln_diamonds_pop" "0.313794535588015" "0.000427916545590119"
[11,] "atlantic_dist_min" "0.725556284996984" "6.49073307389046e-14"
[12,] "indian_dist_min" "0.803117454548449" "3.44740790902477e-16"
[13,] "saharan_dist_min" "0.87224731941177" "6.07174170092653e-19"
[14,] "red_sea_dist_min" "0.871485990898349" "1.06323009420501e-18"
Geary C statistic
[1,] "0.562737299594007" "5.5455869476088e-05"
[2,] "0.551065826441476" "1.47133692708474e-05"
[3,] "0.799389686320364" "0.0686756946949285"
[4,] "0.647335019269922" "0.00217855925159111"
[5,] "0.326450740075291" "2.62229760519281e-10"
[6,] "0.683154733931286" "0.00145453722187371"
[7,] "0.286658983545689" "8.78501766139294e-12"
[8,] "0.801654331091" "0.0291452887728903"
[9,] "0.522042351663106" "8.61703934180597e-06"
[10,] "0.674172952769783" "0.00274798696239173"
[11,] "0.242008049696664" "3.78528596268506e-11"
[12,] "0.205462740496694" "2.38830906084901e-13"
[13,] "0.0968316131722766" "1.04236233334834e-15"
[14,] "0.105249867609864" "2.72806775219181e-16"
showMoran <- function(shf, lmii, tlab, roi)
{
n_colors = 256
my_gradient = c('deepskyblue3', 'snow', 'firebrick3')
palette = colorRampPalette(my_gradient, space = "Lab")(n_colors)
vals = lmii[,1]
#min_v = min(vals); max_v = max(vals);
ref_v = 1;#max(abs(vals))
min_v = -ref_v; max_v = ref_v;
norm_v = (vals-min_v)/(max_v-min_v)
norm_v[norm_v < 0] = 0
norm_v[norm_v > 1] = 1
pcol = palette[round(norm_v * (n_colors-1)) + 1]
#tcol = ifelse(lmii[,5]>0.05, "darkgrey", "black")
tcol = ifelse(1:length(vals) %in% roi, 'black', 'darkgray')
lcex = ifelse(1:length(vals) %in% roi, 0.7, 0.5)
par(mar = c(0.0, 0.0, 1.0, 3))
plot(shf, col=pcol, lwd=0.5, border='olivedrab', main=paste('local Moran: ', tlab))
# Colorbar
image.plot(legend.only = TRUE, zlim = c(min_v, max_v),
col = palette,
axis.args=list(cex.axis=0.8))
coords = coordinates(slaveTrade)
text(coords[,1], coords[,2], labels=ds$isocode, cex=lcex, col=tcol)
par(mar=c(5,4,4,2)+0.1)
}
showClusterMap <- function(ds, lmii, col_lab, significance_level=0.05)
{
par(mar = c(0.0, 0.0, 1.0, 4))
p_value = lmii[,5]
Moran_I = lmii[,1]
llist = c('HH', 'LL', 'HL', 'LH', 'NS')
clist = c('firebrick', 'royalblue', 'indianred1', 'lightblue1', 'lightgray')
lid <- ifelse(p_value < significance_level,
ifelse(Moran_I > 0,
ifelse(ds[[col_lab]] > mean(ds[[col_lab]]), 1, 2),
ifelse(ds[[col_lab]] > mean(ds[[col_lab]]), 3, 4)),
5)
cluster = llist[lid]
ccol = clist[lid]
plot(ds, col=ccol, lwd=0.5, border='snow', main=paste('clusters: ', col_lab))
legend("right", legend = llist, fill = clist, bty='n')
par(mar=c(5,4,4,2)+0.1)
}
computeLocalAvg <- function(dsp, variable_name, wlist)
{
local_averages <- numeric(nrow(dsp)) # Initialize a vector to store the results
for (i in 1:nrow(dsp))
{
neighbors_i <- wlist$neighbours[[i]] #nblist[[i]] # Neighbors of area i
weights_i <- wlist$weights[[i]] #wlist[[i]]
if (length(neighbors_i) > 0)
{ # Check if the area has neighbors (for islands or edge cases)
neighbor_values <- dsp@data[neighbors_i, variable_name] # values of the neighbors
#print(neighbor_values)
local_averages[i] <- sum(weights_i * neighbor_values)/sum(weights_i) # weighted average
}
else
{ # Handle cases with no neighbors (e.g., islands):
local_averages[i] <- dsp@data[i, variable_name]
warning(paste("Area", i, "has no neighbors.")) # or print a message
}
}
return(local_averages)
}
showMapByQuintiles <- function(dsp, variable_name, mydata)
{
dsp_sf = st_as_sf(dsp) # Convert to sf object
dsp_sf$mydata <- mydata
brks=quantile(mydata, seq(0, 1, 0.2), na.rm=T);
brks[1] = floor(brks[1]*10)/10;
brks[length(brks)]=ceiling(brks[length(brks)]*10)/10
spplot(as_Spatial(dsp_sf), 'mydata', main=paste('Quintile cuts: ', variable_name), lwd=0.5,
at = brks,
col = 'lightskyblue4',
# col.regions = colorRampPalette(rev(brewer.pal(11, "Spectral")))(256),
col.regions = matlab.like.hot(16),
par.settings = list(axis.line = list(col = NA)))
}
# brks=quantile(slaveTrade$ln_pcgdp, seq(0, 1, 0.2), na.rm=T);
# brks[1] = floor(brks[1]*10)/10;
# brks[length(brks)]=ceiling(brks[length(brks)]*10)/10
# spplot(slaveTrade, c('ln_pcgdp'), main='ln_pcgdp', lwd=0.5,
# at = brks,
# #col.regions = colorRampPalette(rev(brewer.pal(11, "Spectral")))(256),
# col.regions = matlab.like.hot(16),
# par.settings = list(axis.line = list(col = NA)))
ln_pcgdp_avg = computeLocalAvg(slaveTrade, 'ln_pcgdp', slaveTrade.lw)
dsp_dummy=slaveTrade; dsp_dummy$ln_pcgdp_avg = ln_pcgdp_avg
showData(dsp_dummy,'ln_pcgdp_avg')
showMapByQuintiles(slaveTrade, 'ln_pcgdp_avg', ln_pcgdp_avg)
moran.test(slaveTrade$ln_pcgdp, slaveTrade.lw)
Moran I test under randomisation
data: slaveTrade$ln_pcgdp
weights: slaveTrade.lw
Moran I statistic standard deviate = 4.3832, p-value = 5.847e-06
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.42321222 -0.02222222 0.01032717
corD<-correlog(coordinates(slaveTrade), slaveTrade$ln_pcgdp, method="Moran")
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
corD = cbind(corD, corD[,1:2]); colnames(corD)[5:6] = c('my_Moran', 'my_pvalue')
for (k in 1:nrow(corD))
{
if (k<=length(slaveTrade.nb.lag))
{
nb_lag = slaveTrade.nb.lag[[k]]
lw_klag <- nb2listw(nb_lag, style = 'W')
resk = moran.test(slaveTrade$ln_pcgdp, lw_klag)
corD[k,5] = resk$estimate[1];
corD[k,6] = resk$p.value;
}
else
{
corD[k,5] = NA
corD[k,6] = NA
}
}
(corD)
dist.class coef p.value n my_Moran my_pvalue
[1,] 472379 0.581321114 1.430501e-05 102 0.42321222 5.846936e-06
[2,] 1103686 0.214094436 8.420569e-03 214 -0.09532285 8.516463e-01
[3,] 1734993 -0.145930298 8.950403e-01 218 -0.20962476 9.985113e-01
[4,] 2366300 -0.244325064 9.950600e-01 264 -0.11027096 9.035166e-01
[5,] 2997607 -0.122315364 8.783822e-01 262 0.13563682 2.773817e-02
[6,] 3628914 0.155421053 1.812701e-02 256 NA NA
[7,] 4260221 -0.075661501 6.954287e-01 188 NA NA
[8,] 4891528 -0.285111943 9.964635e-01 208 NA NA
[9,] 5522835 -0.076877245 6.657695e-01 156 NA NA
[10,] 6154142 0.130179157 9.734806e-02 126 NA NA
[11,] 6785449 0.183756900 8.884203e-02 62 NA NA
[12,] 7416756 -0.009985953 2.489103e-01 14 NA NA
plot(x=1:nrow(corD), y=corD[,'coef'], col='black', xlim=c(0,nrow(corD)), bty='n'); grid()
points(x=1:nrow(corD), y=corD[,'my_Moran'], col='firebrick', pch=16, add=T)
Warning: "add" is not a graphical parameter
abline(a=0,b=0,lwd=1)
corD<-correlog(coordinates(slaveTrade), ln_pcgdp_avg, method="Moran")
Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs
par(mar=c(4, 3, 2, 0))
barplot(corD[,'coef'], names.arg=1:length(corD[,'coef']), yaxt = "n",
main=paste('correlog. ', 'ln_pcgdp_avg'), xlab='lag', ylim=c(-1,1),
col='lightskyblue3', border='lightskyblue4')
axis(2, # 2 indicates the left axis (y-axis)
at = seq(-1,1,0.2), # Positions of the ticks
labels = seq(-1,1,0.2))
#grid(lty='solid', lwd=0.5, col='white', nx=NA, ny=NULL)
abline(h = seq(-1,1,0.2), col = "white", lty = 'solid', lwd=0.5) # Add gridlines at y_ticks
par(mar=c(5, 4, 4, 2)+0.1)
NA
NA
set.seed(42)
moran.plot(slaveTrade$ln_pcgdp, slaveTrade.lw,
xlab='ln_pcgdp', asp=1, bty='n')
print(slaveTrade$isocode[c(5,19,23)])
[1] "BWA" "GNQ" "LSO"
lmii = localmoran(slaveTrade$ln_pcgdp, slaveTrade.lw)
# lmii = localmoran_perm(slaveTrade$ln_pcgdp, slaveTrade.lw, nsim=999, iseed=123456789)
showMoran(slaveTrade, lmii, 'ln_pcgdp', c(5,19,23))
showClusterMap(slaveTrade, lmii, 'ln_pcgdp', significance_level=0.05)
ln_export_area_avg = computeLocalAvg(slaveTrade, 'ln_export_area', slaveTrade.lw)
dsp_dummy=slaveTrade; dsp_dummy$ln_export_area_avg = ln_export_area_avg
showData(dsp_dummy,'ln_export_area_avg')
showMapByQuintiles(slaveTrade, 'ln_export_area', ln_export_area_avg)
moran.test(slaveTrade$ln_export_area, slaveTrade.lw)
Moran I test under randomisation
data: slaveTrade$ln_export_area
weights: slaveTrade.lw
Moran I statistic standard deviate = 3.1916, p-value = 0.0007075
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.30551706 -0.02222222 0.01054508
mr = moran.plot(slaveTrade$ln_export_area, slaveTrade.lw,
xlab='ln_export_area', asp=1, bty='n')
print(slaveTrade$isocode[c(40,43)])
[1] "TUN" "ZAF"
lmii = localmoran(slaveTrade$ln_export_area, slaveTrade.lw)
# lmii = localmoran_perm(slaveTrade$ln_export_area, slaveTrade.lw, nsim=999, iseed=123456789)
showMoran(slaveTrade, lmii, 'ln_export_area', c(40,43))
showClusterMap(slaveTrade, lmii, 'ln_export_area')
fmla <- ln_pcgdp ~ ln_export_area + colony1 + colony2 +colony3 + colony4 + colony5 + colony6 + colony7 + rain_min + humid_max + low_temp + ln_coast_area + islam + legor_fr + region_n + ln_gold_pop + ln_oil_pop + ln_diamonds_pop
modLM <- lm(fmla, data = slaveTrade)
summary(modLM)
Call:
lm(formula = fmla, data = slaveTrade)
Residuals:
Min 1Q Median 3Q Max
-0.54419 -0.22749 -0.01176 0.17298 0.75527
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.409676 0.880992 7.276 7.96e-08 ***
ln_export_area -0.090781 0.025850 -3.512 0.00158 **
colony1 0.947924 0.449902 2.107 0.04455 *
colony2 1.043905 0.441921 2.362 0.02563 *
colony3 0.867653 0.454028 1.911 0.06667 .
colony4 0.257498 0.582354 0.442 0.66189
colony5 1.559744 0.705872 2.210 0.03579 *
colony6 1.570997 0.727731 2.159 0.03992 *
colony7 0.469763 0.648013 0.725 0.47473
rain_min 0.001658 0.007255 0.229 0.82096
humid_max 0.016627 0.008935 1.861 0.07367 .
low_temp -0.053279 0.016356 -3.257 0.00303 **
ln_coast_area 0.099378 0.035313 2.814 0.00901 **
islam -0.002329 0.002782 -0.837 0.40971
legor_fr 0.010320 0.421089 0.025 0.98063
region_n -0.131477 0.379578 -0.346 0.73174
ln_gold_pop 0.012940 0.013661 0.947 0.35193
ln_oil_pop 0.077399 0.022307 3.470 0.00177 **
ln_diamonds_pop -0.021472 0.035204 -0.610 0.54701
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.404 on 27 degrees of freedom
Multiple R-squared: 0.8326, Adjusted R-squared: 0.7209
F-statistic: 7.459 on 18 and 27 DF, p-value: 2.204e-06
moran.test(residuals(modLM), slaveTrade.lw)
Moran I test under randomisation
data: residuals(modLM)
weights: slaveTrade.lw
Moran I statistic standard deviate = -0.36075, p-value = 0.6409
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
-0.05897598 -0.02222222 0.01037988
lm.morantest(modLM, slaveTrade.lw)
Global Moran I for regression residuals
data:
model: lm(formula = fmla, data = slaveTrade)
weights: slaveTrade.lw
Moran I statistic standard deviate = 0.55695, p-value = 0.2888
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
-0.058975984 -0.108068846 0.007769563
Anche se dal test di Moran risulta assenza di struttura spaziale sui
residui, e pertanto quanto osservato su ln_pcgdp è già
spiegato dalle covariate impiegate, applichiamo anche il test dei
moltiplicatori di Lagrange (o test di Rao) per valutare la necessità di
un modello a struttura spaziale, tramite lm.RStests().
La funzione lm.RStests() restituisce
#Rao's score (a.k.a Lagrange multiplier) diagnostics
lmtest = lm.RStests(modLM,listw = slaveTrade.lw, test="all" )
summary(lmtest)
Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence
data:
model: lm(formula = fmla, data = slaveTrade)
test weights: slaveTrade.lw
statistic parameter p.value
RSerr 0.3052046 1 0.5806
RSlag 0.0099526 1 0.9205
adjRSerr 0.4906222 1 0.4836
adjRSlag 0.1953702 1 0.6585
SARMA 0.5005748 2 0.7786
Gli elvati p-value forniti da tutti i test sembrano indicare che non esista nessuna dipendenza statisticamente significativa da una struttura spaziale.
Decidiamo comunque di applicare anche gli Spatial Error Model (SEM), Spatial Lag Model (SLM) e lo Spatial Durbin Model (SDM) per ulteriore conferma.
library(spatialreg)
modSEM <- errorsarlm(fmla, data = slaveTrade, listw = slaveTrade.lw)
summary(modSEM)
Call:errorsarlm(formula = fmla, data = slaveTrade, listw = slaveTrade.lw)
Residuals:
Min 1Q Median 3Q Max
-0.492646 -0.216874 -0.047111 0.218542 0.611678
Type: error
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.7163055 0.6755299 6.9816 2.918e-12
ln_export_area -0.0599673 0.0176409 -3.3993 0.0006755
colony1 2.1542058 0.3550998 6.0665 1.307e-09
colony2 1.8236501 0.3169903 5.7530 8.767e-09
colony3 1.2430111 0.3353017 3.7071 0.0002096
colony4 1.1397081 0.3827486 2.9777 0.0029043
colony5 2.5782346 0.4802824 5.3682 7.954e-08
colony6 2.8080166 0.5338434 5.2600 1.441e-07
colony7 1.6344453 0.5735883 2.8495 0.0043787
rain_min 0.0082390 0.0047167 1.7468 0.0806799
humid_max 0.0230702 0.0067576 3.4140 0.0006402
low_temp -0.0772995 0.0108697 -7.1115 1.148e-12
ln_coast_area 0.1478062 0.0235635 6.2727 3.549e-10
islam -0.0022717 0.0017485 -1.2992 0.1938658
legor_fr 0.6507630 0.3589125 1.8132 0.0698084
region_n -0.2911240 0.2332548 -1.2481 0.2119964
ln_gold_pop 0.0158570 0.0106232 1.4927 0.1355246
ln_oil_pop 0.0514046 0.0142822 3.5992 0.0003192
ln_diamonds_pop 0.0016549 0.0251344 0.0658 0.9475043
Lambda: -0.91746, LR test value: 2.344, p-value: 0.12577
Asymptotic standard error: 0.1741
z-value: -5.2696, p-value: 1.3672e-07
Wald statistic: 27.769, p-value: 1.3672e-07
Log likelihood: -10.15064 for error model
ML residual variance (sigma squared): 0.075562, (sigma: 0.27489)
Number of observations: 46
Number of parameters estimated: 21
AIC: 62.301, (AIC for lm: 62.645)
modSEM$coefficients["lambda"]
<NA>
NA
moran.test(residuals(modSEM), slaveTrade.lw)
Moran I test under randomisation
data: residuals(modSEM)
weights: slaveTrade.lw
Moran I statistic standard deviate = -0.27332, p-value = 0.6077
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
-0.05021622 -0.02222222 0.01049021
modSLM <- lagsarlm(fmla, data=slaveTrade, listw=slaveTrade.lw)
summary(modSLM)
Call:lagsarlm(formula = fmla, data = slaveTrade, listw = slaveTrade.lw)
Residuals:
Min 1Q Median 3Q Max
-0.5410864 -0.2299147 -0.0046404 0.1681190 0.7515161
Type: lag
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.5317530 1.4034319 4.6541 3.254e-06
ln_export_area -0.0912573 0.0213286 -4.2786 1.880e-05
colony1 0.9563989 0.3455217 2.7680 0.0056404
colony2 1.0558252 0.3397165 3.1080 0.0018838
colony3 0.8774812 0.3498822 2.5079 0.0121439
colony4 0.2609079 0.4499668 0.5798 0.5620239
colony5 1.5773265 0.5431796 2.9039 0.0036857
colony6 1.5721306 0.5627857 2.7935 0.0052144
colony7 0.4754505 0.4972550 0.9562 0.3389963
rain_min 0.0016561 0.0056141 0.2950 0.7680055
humid_max 0.0166923 0.0068453 2.4385 0.0147478
low_temp -0.0542420 0.0135019 -4.0174 5.885e-05
ln_coast_area 0.1000657 0.0270624 3.6976 0.0002177
islam -0.0023246 0.0021323 -1.0902 0.2756382
legor_fr 0.0114815 0.3240160 0.0354 0.9717329
region_n -0.1394538 0.2951216 -0.4725 0.6365484
ln_gold_pop 0.0129179 0.0104645 1.2345 0.2170335
ln_oil_pop 0.0780160 0.0189756 4.1114 3.933e-05
ln_diamonds_pop -0.0208386 0.0272283 -0.7653 0.4440767
Rho: -0.016867, LR test value: 0.011338, p-value: 0.9152
Asymptotic standard error: 0.14846
z-value: -0.11361, p-value: 0.90955
Wald statistic: 0.012908, p-value: 0.90955
Log likelihood: -11.31695 for lag model
ML residual variance (sigma squared): 0.095761, (sigma: 0.30945)
Number of observations: 46
Number of parameters estimated: 21
AIC: 64.634, (AIC for lm: 62.645)
LM test for residual autocorrelation
test value: 0.55075, p-value: 0.45801
moran.test(residuals(modSLM), slaveTrade.lw)
Moran I test under randomisation
data: residuals(modSLM)
weights: slaveTrade.lw
Moran I statistic standard deviate = -0.31678, p-value = 0.6243
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
-0.05450516 -0.02222222 0.01038525
modSDM <- lagsarlm(fmla, data=slaveTrade, listw=slaveTrade.lw, type='mixed')
summary(modSDM)
Call:lagsarlm(formula = fmla, data = slaveTrade, listw = slaveTrade.lw,
type = "mixed")
Residuals:
Min 1Q Median 3Q Max
-0.3454587 -0.1044586 -0.0089223 0.0780763 0.2933857
Type: mixed
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.02561036 3.09324858 1.3014 0.1931153
ln_export_area -0.15226732 0.02643228 -5.7607 8.379e-09
colony1 1.64439400 0.42314240 3.8861 0.0001018
colony2 1.74069948 0.35216696 4.9428 7.700e-07
colony3 1.03617073 0.36543845 2.8354 0.0045766
colony4 0.20900672 0.43540879 0.4800 0.6312102
colony5 1.65660160 0.38651872 4.2860 1.820e-05
colony6 2.51255171 0.53000474 4.7406 2.131e-06
colony7 1.60481168 0.60658073 2.6457 0.0081530
rain_min -0.02337578 0.00568558 -4.1114 3.932e-05
humid_max 0.03772204 0.00820499 4.5975 4.277e-06
low_temp -0.02299780 0.01526461 -1.5066 0.1319108
ln_coast_area 0.15205182 0.02752320 5.5245 3.304e-08
islam -0.00514750 0.00189003 -2.7235 0.0064594
legor_fr -0.06708916 0.45885552 -0.1462 0.8837558
region_n -0.76731534 0.24208578 -3.1696 0.0015265
ln_gold_pop 0.04240323 0.01125507 3.7675 0.0001649
ln_oil_pop 0.13613948 0.02163123 6.2937 3.101e-10
ln_diamonds_pop -0.06507716 0.02402345 -2.7089 0.0067506
lag.ln_export_area 0.05496841 0.04638199 1.1851 0.2359684
lag.colony1 5.14461052 1.10174121 4.6695 3.019e-06
lag.colony2 4.22291041 0.93596123 4.5118 6.427e-06
lag.colony3 2.10441422 0.91541635 2.2989 0.0215129
lag.colony4 3.44869799 0.85956982 4.0121 6.018e-05
lag.colony5 6.41462091 1.37376814 4.6694 3.021e-06
lag.colony6 1.53701918 1.58396026 0.9704 0.3318647
lag.colony7 7.93422220 1.50120046 5.2853 1.255e-07
lag.rain_min 0.00019509 0.01326682 0.0147 0.9882675
lag.humid_max -0.00484986 0.02456831 -0.1974 0.8435121
lag.low_temp -0.06941803 0.03564798 -1.9473 0.0514964
lag.ln_coast_area 0.12894507 0.06492060 1.9862 0.0470115
lag.islam -0.00049062 0.00427151 -0.1149 0.9085565
lag.legor_fr 1.04046912 1.22216067 0.8513 0.3945828
lag.region_n 0.00018129 0.70397041 0.0003 0.9997945
lag.ln_gold_pop 0.10107588 0.02706225 3.7349 0.0001878
lag.ln_oil_pop -0.16916945 0.04853784 -3.4853 0.0004916
lag.ln_diamonds_pop 0.01798481 0.06894153 0.2609 0.7941924
Rho: -0.48882, LR test value: 3.8581, p-value: 0.049506
Asymptotic standard error: 0.20297
z-value: -2.4084, p-value: 0.016023
Wald statistic: 5.8004, p-value: 0.016023
Log likelihood: 20.26546 for mixed model
ML residual variance (sigma squared): 0.02305, (sigma: 0.15182)
Number of observations: 46
Number of parameters estimated: 39
AIC: 37.469, (AIC for lm: 39.327)
LM test for residual autocorrelation
test value: 11.48, p-value: 0.0007035
moran.test(residuals(modSDM), slaveTrade.lw)
Moran I test under randomisation
data: residuals(modSDM)
weights: slaveTrade.lw
Moran I statistic standard deviate = -0.81573, p-value = 0.7927
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
-0.10527475 -0.02222222 0.01036613
cmpMat = c()
cmpMat[1] = AIC(modLM)
cmpMat[2] = AIC(modSEM)
cmpMat[3] = AIC(modSLM)
cmpMat[4] = AIC(modSDM)
names(cmpMat) = c('LM', 'SEM', 'SLM', 'SDM')
(cmpMat)
LM SEM SLM SDM
62.64525 62.30128 64.63391 37.46908